home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / C / Applications / POV-Ray 3.0.2 / src / SOURCE / POLY.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-04-25  |  27.9 KB  |  1,570 lines  |  [TEXT/CWIE]
  1. /****************************************************************************
  2. *                poly.c
  3. *
  4. *  This module implements the code for general 3 variable polynomial shapes
  5. *
  6. *  This file was written by Alexander Enzmann.  He wrote the code for
  7. *  4th - 6th order shapes and generously provided us these enhancements.
  8. *
  9. *  from Persistence of Vision(tm) Ray Tracer
  10. *  Copyright 1996 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  NOTICE: This source code file is provided so that users may experiment
  13. *  with enhancements to POV-Ray and to port the software to platforms other 
  14. *  than those supported by the POV-Ray Team.  There are strict rules under
  15. *  which you are permitted to use this file.  The rules are in the file
  16. *  named POVLEGAL.DOC which should be distributed with this file. If 
  17. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  19. *  Forum.  The latest version of POV-Ray may be found there as well.
  20. *
  21. * This program is based on the popular DKB raytracer version 2.12.
  22. * DKBTrace was originally written by David K. Buck.
  23. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  24. *
  25. *****************************************************************************/
  26.  
  27. #include "frame.h"
  28. #include "vector.h"
  29. #include "povproto.h"
  30. #include "bbox.h"
  31. #include "polysolv.h"
  32. #include "matrices.h"
  33. #include "objects.h"
  34. #include "poly.h"
  35. #include "povray.h"
  36.  
  37. /*
  38.  * Basic form of a quartic equation:
  39.  *
  40.  *  a00*x^4    + a01*x^3*y   + a02*x^3*z   + a03*x^3     + a04*x^2*y^2 +
  41.  *  a05*x^2*y*z+ a06*x^2*y   + a07*x^2*z^2 + a08*x^2*z   + a09*x^2     +
  42.  *  a10*x*y^3  + a11*x*y^2*z + a12*x*y^2   + a13*x*y*z^2 + a14*x*y*z   +
  43.  *  a15*x*y    + a16*x*z^3   + a17*x*z^2   + a18*x*z     + a19*x       + a20*y^4   +
  44.  *  a21*y^3*z  + a22*y^3     + a23*y^2*z^2 + a24*y^2*z   + a25*y^2     + a26*y*z^3 +
  45.  *  a27*y*z^2  + a28*y*z     + a29*y       + a30*z^4     + a31*z^3     + a32*z^2   + a33*z + a34
  46.  *
  47.  */
  48.  
  49.  
  50.  
  51. /*****************************************************************************
  52. * Local preprocessor defines
  53. ******************************************************************************/
  54.  
  55. #define DEPTH_TOLERANCE 1.0e-4
  56. #define INSIDE_TOLERANCE 1.0e-4
  57. #define ROOT_TOLERANCE 1.0e-4
  58. #define COEFF_LIMIT 1.0e-20
  59. #define BINOMSIZE 40
  60.  
  61.  
  62.  
  63. /*****************************************************************************
  64. * Local typedefs
  65. ******************************************************************************/
  66.  
  67.  
  68.  
  69. /*****************************************************************************
  70. * Static functions
  71. ******************************************************************************/
  72.  
  73. static int intersect PARAMS((RAY *Ray, int Order, DBL *Coeffs, int Sturm_Flag,
  74. DBL *Depths));
  75. static void normal0 PARAMS((VECTOR Result, int Order, DBL *Coeffs,
  76. VECTOR IPoint));
  77. static void normal1 PARAMS((VECTOR Result, int Order, DBL *Coeffs,
  78. VECTOR IPoint));
  79. static DBL inside PARAMS((VECTOR IPoint, int Order, DBL *Coeffs));
  80. static int intersect_linear PARAMS((RAY *ray, DBL *Coeffs, DBL *Depths));
  81. static int intersect_quadratic PARAMS((RAY *ray, DBL *Coeffs, DBL *Depths));
  82. static int factor_out PARAMS((int n, int i, int *c, int *s));
  83. static long binomial PARAMS((int n, int r));
  84. static void factor1 PARAMS((int n, int *c, int *s));
  85.  
  86. /* unused
  87. static DBL evaluate_linear PARAMS((VECTOR P, DBL *a));
  88. static DBL evaluate_quadratic PARAMS((VECTOR P, DBL *a));
  89. */
  90.  
  91. static int All_Poly_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  92. static int Inside_Poly PARAMS((VECTOR IPoint, OBJECT *Object));
  93. static void Poly_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  94. static void *Copy_Poly PARAMS((OBJECT *Object));
  95. static void Translate_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  96. static void Rotate_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  97. static void Scale_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  98. static void Transform_Poly PARAMS((OBJECT *Object, TRANSFORM *Trans));
  99. static void Invert_Poly PARAMS((OBJECT *Object));
  100. static void Destroy_Poly PARAMS((OBJECT *Object));
  101.  
  102. /*****************************************************************************
  103. * Local variables
  104. ******************************************************************************/
  105.  
  106. METHODS Poly_Methods =
  107. {
  108.   All_Poly_Intersections,
  109.   Inside_Poly, Poly_Normal, Copy_Poly,
  110.   Translate_Poly, Rotate_Poly,
  111.   Scale_Poly, Transform_Poly, Invert_Poly, Destroy_Poly
  112. };
  113.  
  114.  
  115.  
  116. /* The following table contains the binomial coefficients up to 15 */
  117.  
  118. static int binomials[15][15] =
  119. {
  120.   {1,  0,  0,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  121.   {1,  1,  0,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  122.   {1,  2,  1,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  123.   {1,  3,  3,  1,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  124.   {1,  4,  6,  4,   1,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  125.   {1,  5, 10, 10,   5,   1,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  126.   {1,  6, 15, 20,  15,   6,   1,   0,   0,   0,   0,  0,  0,  0,  0},
  127.   {1,  7, 21, 35,  35,  21,   7,   1,   0,   0,   0,  0,  0,  0,  0},
  128.   {1,  8, 28, 56,  70,  56,  28,   8,   1,   0,   0,  0,  0,  0,  0},
  129.   {1,  9, 36, 84, 126, 126,  84,  36,   9,   1,   0,  0,  0,  0,  0},
  130.   {1, 10, 45,120, 210, 252, 210, 120,  45,  10,   1,  0,  0,  0,  0},
  131.   {1, 11, 55,165, 330, 462, 462, 330, 165,  55,  11,  1,  0,  0,  0},
  132.   {1, 12, 66,220, 495, 792, 924, 792, 495, 220,  66, 12,  1,  0,  0},
  133.   {1, 13, 78,286, 715,1287,1716,1716,1287, 715, 286, 78, 13,  1,  0},
  134.   {1, 14, 91,364,1001,2002,3003,3432,3003,2002,1001,364, 91, 14,  1}
  135. };
  136.  
  137. static DBL eqn_v[3][MAX_ORDER+1], eqn_vt[3][MAX_ORDER+1];
  138.  
  139.  
  140.  
  141.  
  142. /*****************************************************************************
  143. *
  144. * FUNCTION
  145. *
  146. *   All_Poly_Intersections
  147. *
  148. * INPUT
  149. *   
  150. * OUTPUT
  151. *   
  152. * RETURNS
  153. *   
  154. * AUTHOR
  155. *
  156. *   Alexander Enzmann
  157. *   
  158. * DESCRIPTION
  159. *
  160. *   -
  161. *
  162. * CHANGES
  163. *
  164. *   -
  165. *
  166. ******************************************************************************/
  167.  
  168. static int All_Poly_Intersections(Object, Ray, Depth_Stack)
  169. OBJECT *Object;
  170. RAY *Ray;
  171. ISTACK *Depth_Stack;
  172. {
  173.   POLY *Poly = (POLY *) Object;
  174.   DBL Depths[MAX_ORDER], len;
  175.   VECTOR  IPoint;
  176.   int cnt, i, j, Intersection_Found, same_root;
  177.   RAY New_Ray;
  178.  
  179.   /* Transform the ray into the polynomial's space */
  180.  
  181.   MInvTransPoint(New_Ray.Initial, Ray->Initial, Poly->Trans);
  182.   MInvTransDirection(New_Ray.Direction, Ray->Direction, Poly->Trans);
  183.  
  184.   VLength(len, New_Ray.Direction);
  185.   VInverseScaleEq(New_Ray.Direction, len);
  186.  
  187.   Intersection_Found = FALSE;
  188.  
  189.   Increase_Counter(stats[Ray_Poly_Tests]);
  190.  
  191.   switch (Poly->Order)
  192.   {
  193.     case 1:
  194.  
  195.       cnt = intersect_linear(&New_Ray, Poly->Coeffs, Depths);
  196.  
  197.       break;
  198.  
  199.     case 2:
  200.  
  201.       cnt = intersect_quadratic(&New_Ray, Poly->Coeffs, Depths);
  202.  
  203.       break;
  204.  
  205.     default:
  206.  
  207.       cnt = intersect(&New_Ray, Poly->Order, Poly->Coeffs, Test_Flag(Poly, STURM_FLAG), Depths);
  208.   }
  209.  
  210.   if (cnt > 0)
  211.   {
  212.     Increase_Counter(stats[Ray_Poly_Tests_Succeeded]);
  213.   }
  214.  
  215.   for (i = 0; i < cnt; i++)
  216.   {
  217.     if (Depths[i] > DEPTH_TOLERANCE)
  218.     {
  219.       same_root = FALSE;
  220.  
  221.       for (j = 0; j < i; j++)
  222.       {
  223.         if (Depths[i] == Depths[j])
  224.         {
  225.           same_root = TRUE;
  226.  
  227.           break;
  228.         }
  229.       }
  230.  
  231.       if (!same_root)
  232.       {
  233.         VEvaluateRay(IPoint, New_Ray.Initial, Depths[i], New_Ray.Direction);
  234.  
  235.         /* Transform the point into world space */
  236.  
  237.         MTransPoint(IPoint, IPoint, Poly->Trans);
  238.  
  239.         if (Point_In_Clip(IPoint, Object->Clip))
  240.         {
  241.           push_entry(Depths[i] / len,IPoint,Object,Depth_Stack);
  242.  
  243.           Intersection_Found = TRUE;
  244.         }
  245.       }
  246.     }
  247.   }
  248.  
  249.   return (Intersection_Found);
  250. }
  251.  
  252.  
  253.  
  254. /*****************************************************************************
  255. *
  256. * FUNCTION
  257. *
  258. *   evaluate_linear
  259. *
  260. * INPUT
  261. *   
  262. * OUTPUT
  263. *   
  264. * RETURNS
  265. *   
  266. * AUTHOR
  267. *
  268. *   Alexander Enzmann
  269. *   
  270. * DESCRIPTION
  271. *
  272. *   -
  273. *
  274. * CHANGES
  275. *
  276. *   -
  277. *
  278. ******************************************************************************/
  279.  
  280. /* For speedup of low order polynomials, expand out the terms
  281.    involved in evaluating the poly. */
  282. /* unused
  283. static DBL evaluate_linear(P, a)
  284. VECTOR P;
  285. DBL *a;
  286. {
  287.   return(a[0] * P[X]) + (a[1] * P[Y]) + (a[2] * P[Z]) + a[3];
  288. }
  289. */
  290.  
  291.  
  292.  
  293. /*****************************************************************************
  294. *
  295. * FUNCTION
  296. *
  297. *   evaluate_quadratic
  298. *
  299. * INPUT
  300. *   
  301. * OUTPUT
  302. *   
  303. * RETURNS
  304. *   
  305. * AUTHOR
  306. *
  307. *   Alexander Enzmann
  308. *   
  309. * DESCRIPTION
  310. *
  311. *   -
  312. *
  313. * CHANGES
  314. *
  315. *   -
  316. *
  317. ******************************************************************************/
  318.  
  319. /*
  320. static DBL evaluate_quadratic(P, a)
  321. VECTOR P;
  322. DBL *a;
  323. {
  324.   DBL x, y, z;
  325.  
  326.   x = P[X];
  327.   y = P[Y];
  328.   z = P[Z];
  329.  
  330.   return(a[0] * x * x + a[1] * x * y + a[2] * x * z +
  331.          a[3] * x     + a[4] * y * y + a[5] * y * z +
  332.          a[6] * y     + a[7] * z * z + a[8] * z     + a[9]);
  333. }
  334. */
  335.  
  336.  
  337.  
  338. /*****************************************************************************
  339. *
  340. * FUNCTION
  341. *
  342. *   factor_out
  343. *
  344. * INPUT
  345. *   
  346. * OUTPUT
  347. *   
  348. * RETURNS
  349. *   
  350. * AUTHOR
  351. *
  352. *   Alexander Enzmann
  353. *   
  354. * DESCRIPTION
  355. *
  356. *   Remove all factors of i from n.
  357. *
  358. * CHANGES
  359. *
  360. *   -
  361. *
  362. ******************************************************************************/
  363.  
  364. static int factor_out(n, i, c, s)
  365. int n, i, *c, *s;
  366. {
  367.   while (!(n % i))
  368.   {
  369.     n /= i;
  370.  
  371.     s[(*c)++] = i;
  372.   }
  373.  
  374.   return(n);
  375. }
  376.  
  377.  
  378.  
  379. /*****************************************************************************
  380. *
  381. * FUNCTION
  382. *
  383. *   factor1
  384. *
  385. * INPUT
  386. *   
  387. * OUTPUT
  388. *   
  389. * RETURNS
  390. *   
  391. * AUTHOR
  392. *
  393. *   Alexander Enzmann
  394. *   
  395. * DESCRIPTION
  396. *
  397. *   Find all prime factors of n. (Note that n must be less than 2^15.
  398. *
  399. * CHANGES
  400. *
  401. *   -
  402. *
  403. ******************************************************************************/
  404.  
  405. static void factor1(n, c, s)
  406. int n, *c, *s;
  407. {
  408.   int i,k;
  409.  
  410.   /* First factor out any 2s. */
  411.  
  412.   n = factor_out(n, 2, c, s);
  413.  
  414.   /* Now any odd factors. */
  415.  
  416.   k = (int)sqrt((DBL)n) + 1;
  417.  
  418.   for (i = 3; (n > 1) && (i <= k); i += 2)
  419.   {
  420.     if (!(n % i))
  421.     {
  422.       n = factor_out(n, i, c, s);
  423.       k = (int)sqrt((DBL)n)+1;
  424.     }
  425.   }
  426.  
  427.   if (n > 1)
  428.   {
  429.     s[(*c)++] = n;
  430.   }
  431. }
  432.  
  433.  
  434.  
  435. /*****************************************************************************
  436. *
  437. * FUNCTION
  438. *
  439. *   binomial
  440. *
  441. * INPUT
  442. *   
  443. * OUTPUT
  444. *   
  445. * RETURNS
  446. *   
  447. * AUTHOR
  448. *
  449. *   Alexander Enzmann
  450. *   
  451. * DESCRIPTION
  452. *
  453. *   Calculate the binomial coefficent of n, r.
  454. *
  455. * CHANGES
  456. *
  457. *   -
  458. *
  459. ******************************************************************************/
  460.  
  461. static long binomial(n, r)
  462. int n, r;
  463. {
  464.   int h,i,j,k,l;
  465.   unsigned long result;
  466.   static int stack1[BINOMSIZE], stack2[BINOMSIZE];
  467.  
  468.   if ((n < 0) || (r < 0) || (r > n))
  469.   {
  470.     result = 0L;
  471.   }
  472.   else
  473.   {
  474.     if (r == n)
  475.     {
  476.       result = 1L;
  477.     }
  478.     else
  479.     {
  480.       if ((r < 15) && (n < 15))
  481.       {
  482.         result = (long)binomials[n][r];
  483.       }
  484.       else
  485.       {
  486.         j = 0;
  487.  
  488.         for (i = r + 1; i <= n; i++)
  489.         {
  490.           stack1[j++] = i;
  491.         }
  492.  
  493.         for (i = 2; i <= (n-r); i++)
  494.         {
  495.           h = 0;
  496.  
  497.           factor1(i, &h, stack2);
  498.  
  499.           for (k = 0; k < h; k++)
  500.           {
  501.             for (l = 0; l < j; l++)
  502.             {
  503.               if (!(stack1[l] % stack2[k]))
  504.               {
  505.                 stack1[l] /= stack2[k];
  506.  
  507.                 goto l1;
  508.               }
  509.             }
  510.           }
  511.  
  512.           /* Error if we get here */
  513. /*        Debug_Info("Failed to factor\n");*/
  514. l1:;
  515.         }
  516.  
  517.         result = 1;
  518.  
  519.         for (i = 0; i < j; i++)
  520.         {
  521.           result *= stack1[i];
  522.         }
  523.       }
  524.     }
  525.   }
  526.  
  527.   return(result);
  528. }
  529.  
  530.  
  531.  
  532. /*****************************************************************************
  533. *
  534. * FUNCTION
  535. *
  536. *   inside
  537. *
  538. * INPUT
  539. *   
  540. * OUTPUT
  541. *   
  542. * RETURNS
  543. *   
  544. * AUTHOR
  545. *
  546. *   Alexander Enzmann
  547. *   
  548. * DESCRIPTION
  549. *
  550. *   -
  551. *
  552. * CHANGES
  553. *
  554. *   -
  555. *
  556. ******************************************************************************/
  557.  
  558. static DBL inside(IPoint, Order, Coeffs)
  559. VECTOR  IPoint;
  560. int Order;
  561. DBL *Coeffs;
  562. {
  563.   DBL x[MAX_ORDER+1], y[MAX_ORDER+1];
  564.   DBL z[MAX_ORDER+1], c, Result;
  565.   int i, j, k, term;
  566.  
  567.   x[0] = 1.0;       y[0] = 1.0;       z[0] = 1.0;
  568.   x[1] = IPoint[X]; y[1] = IPoint[Y]; z[1] = IPoint[Z];
  569.  
  570.   for (i = 2; i <= Order; i++)
  571.   {
  572.     x[i] = x[1] * x[i-1];
  573.     y[i] = y[1] * y[i-1];
  574.     z[i] = z[1] * z[i-1];
  575.   }
  576.  
  577.   Result = 0.0;
  578.  
  579.   term = 0;
  580.  
  581.   for (i = Order; i >= 0; i--)
  582.   {
  583.     for (j=Order-i;j>=0;j--)
  584.     {
  585.       for (k=Order-(i+j);k>=0;k--)
  586.       {
  587.         if ((c = Coeffs[term]) != 0.0)
  588.         {
  589.           Result += c * x[i] * y[j] * z[k];
  590.         }
  591.         term++;
  592.       }
  593.     }
  594.   }
  595.  
  596.   return(Result);
  597. }
  598.  
  599.  
  600.  
  601. /*****************************************************************************
  602. *
  603. * FUNCTION
  604. *
  605. *   intersect
  606. *
  607. * INPUT
  608. *   
  609. * OUTPUT
  610. *   
  611. * RETURNS
  612. *   
  613. * AUTHOR
  614. *
  615. *   Alexander Enzmann
  616. *   
  617. * DESCRIPTION
  618. *
  619. *   Intersection of a ray and an arbitrary polynomial function.
  620. *
  621. * CHANGES
  622. *
  623. *   -
  624. *
  625. ******************************************************************************/
  626.  
  627. static int intersect(ray, Order, Coeffs, Sturm_Flag, Depths)
  628. RAY *ray;
  629. int Order, Sturm_Flag;
  630. DBL *Coeffs, *Depths;
  631. {
  632.   DBL eqn[MAX_ORDER+1];
  633.   DBL t[3][MAX_ORDER+1];
  634.   VECTOR  P, D;
  635.   DBL val;
  636.   int h, i, j, k, i1, j1, k1, term;
  637.   int offset;
  638.  
  639.   /* First we calculate the values of the individual powers
  640.      of x, y, and z as they are represented by the ray */
  641.  
  642.   Assign_Vector(P,ray->Initial);
  643.   Assign_Vector(D,ray->Direction);
  644.  
  645.   for (i = 0; i < 3; i++)
  646.   {
  647.     eqn_v[i][0]  = 1.0;
  648.     eqn_vt[i][0] = 1.0;
  649.   }
  650.  
  651.   eqn_v[0][1] = P[X];
  652.   eqn_v[1][1] = P[Y];
  653.   eqn_v[2][1] = P[Z];
  654.  
  655.   eqn_vt[0][1] = D[X];
  656.   eqn_vt[1][1] = D[Y];
  657.   eqn_vt[2][1] = D[Z];
  658.  
  659.   for (i = 2; i <= Order; i++)
  660.   {
  661.     for (j=0;j<3;j++)
  662.     {
  663.      eqn_v[j][i]  = eqn_v[j][1] * eqn_v[j][i-1];
  664.      eqn_vt[j][i] = eqn_vt[j][1] * eqn_vt[j][i-1];
  665.     }
  666.   }
  667.  
  668.   for (i = 0; i <= Order; i++)
  669.   {
  670.     eqn[i] = 0.0;
  671.   }
  672.  
  673.   /* Now walk through the terms of the polynomial.  As we go
  674.      we substitute the ray equation for each of the variables. */
  675.  
  676.   term = 0;
  677.  
  678.   for (i = Order; i >= 0; i--)
  679.   {
  680.     for (h = 0; h <= i; h++)
  681.     {
  682.       t[0][h] = binomial(i, h) * eqn_vt[0][i-h] * eqn_v[0][h];
  683.     }
  684.  
  685.     for (j = Order-i; j >= 0; j--)
  686.     {
  687.       for (h = 0; h <= j; h++)
  688.       {
  689.         t[1][h] = binomial(j, h) * eqn_vt[1][j-h] * eqn_v[1][h];
  690.       }
  691.  
  692.       for (k = Order-(i+j); k >= 0; k--)
  693.       {
  694.         if (Coeffs[term] != 0)
  695.         {
  696.           for (h = 0; h <= k; h++)
  697.           {
  698.             t[2][h] = binomial(k, h) * eqn_vt[2][k-h] * eqn_v[2][h];
  699.           }
  700.  
  701.           /* Multiply together the three polynomials. */
  702.  
  703.           offset = Order - (i + j + k);
  704.  
  705.           for (i1 = 0; i1 <= i; i1++)
  706.           {
  707.             for (j1=0;j1<=j;j1++)
  708.             {
  709.               for (k1=0;k1<=k;k1++)
  710.               {
  711.                 h = offset + i1 + j1 + k1;
  712.                 val = Coeffs[term];
  713.                 val *= t[0][i1];
  714.                 val *= t[1][j1];
  715.                 val *= t[2][k1];
  716.                 eqn[h] += val;
  717.               }
  718.             }
  719.           }
  720.         }
  721.  
  722.         term++;
  723.       }
  724.     }
  725.   }
  726.  
  727.   for (i = 0, j = Order; i <= Order; i++)
  728.   {
  729.     if (eqn[i] != 0.0)
  730.     {
  731.       break;
  732.     }
  733.     else
  734.     {
  735.       j--;
  736.     }
  737.   }
  738.  
  739.   if (j > 1)
  740.   {
  741.     return(Solve_Polynomial(j, &eqn[i], Depths, Sturm_Flag, ROOT_TOLERANCE));
  742.   }
  743.   else
  744.   {
  745.     return(0);
  746.   }
  747. }
  748.  
  749.  
  750.  
  751. /*****************************************************************************
  752. *
  753. * FUNCTION
  754. *
  755. *   intersect_linear
  756. *
  757. * INPUT
  758. *   
  759. * OUTPUT
  760. *   
  761. * RETURNS
  762. *   
  763. * AUTHOR
  764. *
  765. *   Alexander Enzmann
  766. *   
  767. * DESCRIPTION
  768. *
  769. *   Intersection of a ray and a quadratic.
  770. *
  771. * CHANGES
  772. *
  773. *   -
  774. *
  775. ******************************************************************************/
  776.  
  777. static int intersect_linear(ray, Coeffs, Depths)
  778. RAY *ray;
  779. DBL *Coeffs, *Depths;
  780. {
  781.   DBL t0, t1, *a = Coeffs;
  782.  
  783.   t0 = a[0] * ray->Initial[X] + a[1] * ray->Initial[Y] + a[2] * ray->Initial[Z];
  784.   t1 = a[0] * ray->Direction[X] + a[1] * ray->Direction[Y] +
  785.  
  786.   a[2] * ray->Direction[Z];
  787.  
  788.   if (fabs(t1) < EPSILON)
  789.   {
  790.     return(0);
  791.   }
  792.  
  793.   Depths[0] = -(a[3] + t0) / t1;
  794.  
  795.   return(1);
  796. }
  797.  
  798.  
  799.  
  800. /*****************************************************************************
  801. *
  802. * FUNCTION
  803. *
  804. *   intersect_quadratic
  805. *
  806. * INPUT
  807. *   
  808. * OUTPUT
  809. *   
  810. * RETURNS
  811. *   
  812. * AUTHOR
  813. *
  814. *   Alexander Enzmann
  815. *   
  816. * DESCRIPTION
  817. *
  818. *   Intersection of a ray and a quadratic.
  819. *
  820. * CHANGES
  821. *
  822. *   -
  823. *
  824. ******************************************************************************/
  825.  
  826. static int intersect_quadratic(ray, Coeffs, Depths)
  827. RAY *ray;
  828. DBL *Coeffs, *Depths;
  829. {
  830.   DBL x, y, z, x2, y2, z2;
  831.   DBL xx, yy, zz, xx2, yy2, zz2;
  832.   DBL *a, ac, bc, cc, d, t;
  833.  
  834.   x  = ray->Initial[X];
  835.   y  = ray->Initial[Y];
  836.   z  = ray->Initial[Z];
  837.  
  838.   xx = ray->Direction[X];
  839.   yy = ray->Direction[Y];
  840.   zz = ray->Direction[Z];
  841.  
  842.   x2  = x * x;    y2  = y * y;    z2 = z * z;
  843.   xx2 = xx * xx;  yy2 = yy * yy;  zz2 = zz * zz;
  844.  
  845.   a = Coeffs;
  846.  
  847.   /*
  848.    * Determine the coefficients of t^n, where the line is represented
  849.    * as (x,y,z) + (xx,yy,zz)*t.
  850.    */
  851.  
  852.   ac = (a[0]*xx2 + a[1]*xx*yy + a[2]*xx*zz + a[4]*yy2 + a[5]*yy*zz + a[7]*zz2);
  853.  
  854.   bc = (2*a[0]*x*xx + a[1]*(x*yy + xx*y) + a[2]*(x*zz + xx*z) +
  855.         a[3]*xx + 2*a[4]*y*yy + a[5]*(y*zz + yy*z) + a[6]*yy +
  856.         2*a[7]*z*zz + a[8]*zz);
  857.  
  858.   cc = a[0]*x2 + a[1]*x*y + a[2]*x*z + a[3]*x + a[4]*y2 +
  859.        a[5]*y*z + a[6]*y + a[7]*z2 + a[8]*z + a[9];
  860.  
  861.   if (fabs(ac) < COEFF_LIMIT)
  862.   {
  863.     if (fabs(bc) < COEFF_LIMIT)
  864.     {
  865.       return(0);
  866.     }
  867.  
  868.     Depths[0] = -cc / bc;
  869.  
  870.     return(1);
  871.   }
  872.  
  873.   /*
  874.    * Solve the quadratic formula & return results that are
  875.    * within the correct interval.
  876.    */
  877.  
  878.   d = bc * bc - 4.0 * ac * cc;
  879.  
  880.   if (d < 0.0)
  881.   {
  882.     return(0);
  883.   }
  884.  
  885.   d = sqrt(d);
  886.  
  887.   bc = -bc;
  888.  
  889.   t = 2.0 * ac;
  890.  
  891.   Depths[0] = (bc + d) / t;
  892.   Depths[1] = (bc - d) / t;
  893.  
  894.   return(2);
  895. }
  896.  
  897.  
  898.  
  899. /*****************************************************************************
  900. *
  901. * FUNCTION
  902. *
  903. *   normal0
  904. *
  905. * INPUT
  906. *   
  907. * OUTPUT
  908. *   
  909. * RETURNS
  910. *   
  911. * AUTHOR
  912. *
  913. *   Alexander Enzmann
  914. *   
  915. * DESCRIPTION
  916. *
  917. *   Normal to a polynomial (used for polynomials with order > 4).
  918. *
  919. * CHANGES
  920. *
  921. *   -
  922. *
  923. ******************************************************************************/
  924.  
  925. static void normal0(Result, Order, Coeffs, IPoint)
  926. VECTOR Result;
  927. int Order;
  928. DBL *Coeffs;
  929. VECTOR IPoint;
  930.  {
  931.   int i, j, k, term;
  932.   DBL val, *a, x[MAX_ORDER+1], y[MAX_ORDER+1], z[MAX_ORDER+1];
  933.  
  934.   x[0] = 1.0;
  935.   y[0] = 1.0;
  936.   z[0] = 1.0;
  937.  
  938.   x[1] = IPoint[X];
  939.   y[1] = IPoint[Y];
  940.   z[1] = IPoint[Z];
  941.  
  942.   for (i = 2; i <= Order; i++)
  943.   {
  944.     x[i] = IPoint[X] * x[i-1];
  945.     y[i] = IPoint[Y] * y[i-1];
  946.     z[i] = IPoint[Z] * z[i-1];
  947.   }
  948.  
  949.   a = Coeffs;
  950.  
  951.   term = 0;
  952.  
  953.   Make_Vector(Result, 0.0, 0.0, 0.0);
  954.  
  955.   for (i = Order; i >= 0; i--)
  956.   {
  957.     for (j = Order-i; j >= 0; j--)
  958.     {
  959.       for (k = Order-(i+j); k >= 0; k--)
  960.       {
  961.         if (i >= 1)
  962.         {
  963.           val = x[i-1] * y[j] * z[k];
  964.           Result[X] += i * a[term] * val;
  965.         }
  966.  
  967.         if (j >= 1)
  968.         {
  969.           val = x[i] * y[j-1] * z[k];
  970.           Result[Y] += j * a[term] * val;
  971.         }
  972.  
  973.         if (k >= 1)
  974.         {
  975.           val = x[i] * y[j] * z[k-1];
  976.           Result[Z] += k * a[term] * val;
  977.         }
  978.  
  979.         term++;
  980.       }
  981.     }
  982.   }
  983. }
  984.  
  985.  
  986.  
  987. /*****************************************************************************
  988. *
  989. * FUNCTION
  990. *
  991. *   nromal1
  992. *
  993. * INPUT
  994. *   
  995. * OUTPUT
  996. *   
  997. * RETURNS
  998. *   
  999. * AUTHOR
  1000. *
  1001. *   Alexander Enzmann
  1002. *   
  1003. * DESCRIPTION
  1004. *
  1005. *   Normal to a polynomial (for polynomials of order <= 4).
  1006. *
  1007. * CHANGES
  1008. *
  1009. *   -
  1010. *
  1011. ******************************************************************************/
  1012.  
  1013. static void normal1(Result, Order, Coeffs, IPoint)
  1014. VECTOR Result;
  1015. int Order;
  1016. DBL *Coeffs;
  1017. VECTOR IPoint;
  1018. {
  1019.   DBL *a, x, y, z, x2, y2, z2, x3, y3, z3;
  1020.  
  1021.   a = Coeffs;
  1022.  
  1023.   x = IPoint[X];
  1024.   y = IPoint[Y];
  1025.   z = IPoint[Z];
  1026.  
  1027.   switch (Order)
  1028.   {
  1029.     case 1:
  1030.  
  1031.       /* Linear partial derivatives */
  1032.  
  1033.       Make_Vector(Result, a[0], a[1], a[2])
  1034.  
  1035.       break;
  1036.  
  1037.     case 2:
  1038.  
  1039.       /* Quadratic partial derivatives */
  1040.  
  1041.       Result[X] = 2*a[0]*x+a[1]*y+a[2]*z+a[3];
  1042.       Result[Y] = a[1]*x+2*a[4]*y+a[5]*z+a[6];
  1043.       Result[Z] = a[2]*x+a[5]*y+2*a[7]*z+a[8];
  1044.  
  1045.       break;
  1046.  
  1047.     case 3:
  1048.  
  1049.         x2 = x * x;  y2 = y * y;  z2 = z * z;
  1050.  
  1051.         /* Cubic partial derivatives */
  1052.  
  1053.         Result[X] = 3*a[0]*x2 + 2*x*(a[1]*y + a[2]*z + a[3]) + a[4]*y2 +
  1054.                     y*(a[5]*z + a[6]) + a[7]*z2 + a[8]*z + a[9];
  1055.         Result[Y] = a[1]*x2 + x*(2*a[4]*y + a[5]*z + a[6]) + 3*a[10]*y2 +
  1056.                     2*y*(a[11]*z + a[12]) + a[13]*z2 + a[14]*z + a[15];
  1057.         Result[Z] = a[2]*x2 + x*(a[5]*y + 2*a[7]*z + a[8]) + a[11]*y2 +
  1058.                     y*(2*a[13]*z + a[14]) + 3*a[16]*z2 + 2*a[17]*z + a[18];
  1059.  
  1060.         break;
  1061.  
  1062.     case 4:
  1063.  
  1064.       /* Quartic partial derivatives */
  1065.  
  1066.       x2 = x * x;  y2 = y * y;  z2 = z * z;
  1067.       x3 = x * x2; y3 = y * y2; z3 = z * z2;
  1068.  
  1069.       Result[X] = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+
  1070.                   2*x*(a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9])+
  1071.                   a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+
  1072.                   a[16]*z3+a[17]*z2+a[18]*z+a[19];
  1073.       Result[Y] = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+
  1074.                   x*(3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15])+
  1075.                   4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+
  1076.                   a[26]*z3+a[27]*z2+a[28]*z+a[29];
  1077.       Result[Z] = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+
  1078.                   x*(a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18])+
  1079.                   a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+
  1080.                   4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33];
  1081.   }
  1082. }
  1083.  
  1084.  
  1085.  
  1086. /*****************************************************************************
  1087. *
  1088. * FUNCTION
  1089. *
  1090. *   Inside_Poly
  1091. *
  1092. * INPUT
  1093. *   
  1094. * OUTPUT
  1095. *   
  1096. * RETURNS
  1097. *   
  1098. * AUTHOR
  1099. *
  1100. *   Alexander Enzmann
  1101. *
  1102. * DESCRIPTION
  1103. *
  1104. *   -
  1105. *
  1106. * CHANGES
  1107. *
  1108. *   -
  1109. *
  1110. ******************************************************************************/
  1111.  
  1112. static int Inside_Poly (IPoint, Object)
  1113. VECTOR IPoint;
  1114. OBJECT *Object;
  1115. {
  1116.   VECTOR  New_Point;
  1117.   DBL Result;
  1118.  
  1119.   /* Transform the point into polynomial's space */
  1120.  
  1121.   MInvTransPoint(New_Point, IPoint, ((POLY *)Object)->Trans);
  1122.  
  1123.   Result = inside(New_Point, ((POLY *)Object)->Order, ((POLY *)Object)->Coeffs);
  1124.  
  1125.   if (Result < INSIDE_TOLERANCE)
  1126.   {
  1127.     return ((int)(!Test_Flag(Object, INVERTED_FLAG)));
  1128.   }
  1129.   else
  1130.   {
  1131.     return ((int)Test_Flag(Object, INVERTED_FLAG));
  1132.   }
  1133. }
  1134.  
  1135.  
  1136.  
  1137. /*****************************************************************************
  1138. *
  1139. * FUNCTION
  1140. *
  1141. *   Poly_Normal
  1142. *
  1143. * INPUT
  1144. *   
  1145. * OUTPUT
  1146. *
  1147. * RETURNS
  1148. *   
  1149. * AUTHOR
  1150. *
  1151. *   Alexander Enzmann
  1152. *   
  1153. * DESCRIPTION
  1154. *
  1155. *   Normal to a polynomial.
  1156. *
  1157. * CHANGES
  1158. *
  1159. *   -
  1160. *
  1161. ******************************************************************************/
  1162.  
  1163. static void Poly_Normal(Result, Object, Inter)
  1164. OBJECT *Object;
  1165. VECTOR Result;
  1166. INTERSECTION *Inter;
  1167. {
  1168.   DBL val;
  1169.   VECTOR  New_Point;
  1170.   POLY *Shape = (POLY *)Object;
  1171.  
  1172.   /* Transform the point into the polynomials space. */
  1173.  
  1174.   MInvTransPoint(New_Point, Inter->IPoint, Shape->Trans);
  1175.  
  1176.   if (((POLY *)Object)->Order > 4)
  1177.   {
  1178.     normal0(Result, Shape->Order, Shape->Coeffs, New_Point);
  1179.   }
  1180.   else
  1181.   {
  1182.     normal1(Result, Shape->Order, Shape->Coeffs, New_Point);
  1183.   }
  1184.  
  1185.   /* Transform back to world space. */
  1186.  
  1187.   MTransNormal(Result, Result, Shape->Trans);
  1188.  
  1189.   /* Normalize (accounting for the possibility of a 0 length normal). */
  1190.  
  1191.   VDot(val, Result, Result);
  1192.  
  1193.   if (val > 0.0)
  1194.   {
  1195.     val = 1.0 / sqrt(val);
  1196.  
  1197.     VScaleEq(Result, val);
  1198.   }
  1199.   else
  1200.   {
  1201.     Make_Vector(Result, 1.0, 0.0, 0.0)
  1202.   }
  1203. }
  1204.  
  1205.  
  1206.  
  1207. /*****************************************************************************
  1208. *
  1209. * FUNCTION
  1210. *
  1211. *   Translate_Poly
  1212. *
  1213. * INPUT
  1214. *   
  1215. * OUTPUT
  1216. *   
  1217. * RETURNS
  1218. *   
  1219. * AUTHOR
  1220. *
  1221. *   Alexander Enzmann
  1222. *   
  1223. * DESCRIPTION
  1224. *
  1225. *   -
  1226. *
  1227. * CHANGES
  1228. *
  1229. *   -
  1230. *
  1231. ******************************************************************************/
  1232.  
  1233. static void Translate_Poly (Object, Vector, Trans)
  1234. OBJECT *Object;
  1235. VECTOR Vector;
  1236. TRANSFORM *Trans;
  1237. {
  1238.   Transform_Poly(Object, Trans);
  1239. }
  1240.  
  1241.  
  1242.  
  1243. /*****************************************************************************
  1244. *
  1245. * FUNCTION
  1246. *
  1247. *   Rotate_Poly
  1248. *
  1249. * INPUT
  1250. *   
  1251. * OUTPUT
  1252. *   
  1253. * RETURNS
  1254. *   
  1255. * AUTHOR
  1256. *
  1257. *   Alexander Enzmann
  1258. *   
  1259. * DESCRIPTION
  1260. *
  1261. *   -
  1262. *
  1263. * CHANGES
  1264. *
  1265. *   -
  1266. *
  1267. ******************************************************************************/
  1268.  
  1269. static void Rotate_Poly (Object, Vector, Trans)
  1270. OBJECT *Object;
  1271. VECTOR Vector;
  1272. TRANSFORM *Trans;
  1273. {
  1274.   Transform_Poly(Object, Trans);
  1275. }
  1276.  
  1277.  
  1278.  
  1279. /*****************************************************************************
  1280. *
  1281. * FUNCTION
  1282. *
  1283. *   Scale_Poly
  1284. *
  1285. * INPUT
  1286. *   
  1287. * OUTPUT
  1288. *   
  1289. * RETURNS
  1290. *   
  1291. * AUTHOR
  1292. *
  1293. *   Alexander Enzmann
  1294. *   
  1295. * DESCRIPTION
  1296. *
  1297. *   -
  1298. *
  1299. * CHANGES
  1300. *
  1301. *   -
  1302. *
  1303. ******************************************************************************/
  1304.  
  1305. static void Scale_Poly (Object, Vector, Trans)
  1306. OBJECT *Object;
  1307. VECTOR Vector;
  1308. TRANSFORM *Trans;
  1309. {
  1310.   Transform_Poly(Object, Trans);
  1311. }
  1312.  
  1313.  
  1314.  
  1315. /*****************************************************************************
  1316. *
  1317. * FUNCTION
  1318. *
  1319. *   Transform_Poly
  1320. *
  1321. * INPUT
  1322. *   
  1323. * OUTPUT
  1324. *   
  1325. * RETURNS
  1326. *   
  1327. * AUTHOR
  1328. *
  1329. *   Alexander Enzmann
  1330. *   
  1331. * DESCRIPTION
  1332. *
  1333. *   -
  1334. *
  1335. * CHANGES
  1336. *
  1337. *   -
  1338. *
  1339. ******************************************************************************/
  1340.  
  1341. static void Transform_Poly(Object,Trans)
  1342. OBJECT *Object;
  1343. TRANSFORM *Trans;
  1344. {
  1345.   Compose_Transforms(((POLY *)Object)->Trans, Trans);
  1346.  
  1347.   Compute_Poly_BBox((POLY *)Object);
  1348. }
  1349.  
  1350.  
  1351.  
  1352. /*****************************************************************************
  1353. *
  1354. * FUNCTION
  1355. *
  1356. *   Invert_Poly
  1357. *
  1358. * INPUT
  1359. *   
  1360. * OUTPUT
  1361. *   
  1362. * RETURNS
  1363. *   
  1364. * AUTHOR
  1365. *
  1366. *   Alexander Enzmann
  1367. *   
  1368. * DESCRIPTION
  1369. *
  1370. *   -
  1371. *
  1372. * CHANGES
  1373. *
  1374. *   -
  1375. *
  1376. ******************************************************************************/
  1377.  
  1378. static void Invert_Poly(Object)
  1379. OBJECT *Object;
  1380. {
  1381.   Invert_Flag(Object, INVERTED_FLAG);
  1382. }
  1383.  
  1384.  
  1385.  
  1386. /*****************************************************************************
  1387. *
  1388. * FUNCTION
  1389. *
  1390. *   Create_Poly
  1391. *
  1392. * INPUT
  1393. *   
  1394. * OUTPUT
  1395. *   
  1396. * RETURNS
  1397. *   
  1398. * AUTHOR
  1399. *
  1400. *   Alexander Enzmann
  1401. *   
  1402. * DESCRIPTION
  1403. *
  1404. *   -
  1405. *
  1406. * CHANGES
  1407. *
  1408. *   -
  1409. *
  1410. ******************************************************************************/
  1411.  
  1412. POLY *Create_Poly(Order)
  1413. int Order;
  1414. {
  1415.   POLY *New;
  1416.   int i;
  1417.  
  1418.   New = (POLY *)POV_MALLOC(sizeof (POLY), "poly");
  1419.  
  1420.   INIT_OBJECT_FIELDS(New,POLY_OBJECT, &Poly_Methods);
  1421.  
  1422.   New->Order = Order;
  1423.  
  1424.   New->Trans = Create_Transform();
  1425.  
  1426.   New->Coeffs = (DBL *)POV_MALLOC(term_counts(Order) * sizeof(DBL), "coefficients for POLY");
  1427.  
  1428.   for (i = 0; i < term_counts(Order); i++)
  1429.   {
  1430.     New->Coeffs[i] = 0.0;
  1431.   }
  1432.  
  1433.   return(New);
  1434. }
  1435.  
  1436.  
  1437.  
  1438. /*****************************************************************************
  1439. *
  1440. * FUNCTION
  1441. *
  1442. *   Copy_Poly
  1443. *
  1444. * INPUT
  1445. *   
  1446. * OUTPUT
  1447. *   
  1448. * RETURNS
  1449. *   
  1450. * AUTHOR
  1451. *
  1452. *   Alexander Enzmann
  1453. *   
  1454. * DESCRIPTION
  1455. *
  1456. *   Make a copy of a polynomial object.
  1457. *
  1458. * CHANGES
  1459. *
  1460. *   -
  1461. *
  1462. ******************************************************************************/
  1463.  
  1464. static void *Copy_Poly(Object)
  1465. OBJECT *Object;
  1466. {
  1467.   POLY *New;
  1468.   int i;
  1469.  
  1470.   New = Create_Poly (((POLY *)Object)->Order);
  1471.  
  1472.   /* Get rid of transform created in Create_Poly. */
  1473.  
  1474.   Destroy_Transform(New->Trans);
  1475.  
  1476.   Copy_Flag(New, Object, STURM_FLAG);
  1477.   Copy_Flag(New, Object, INVERTED_FLAG);
  1478.  
  1479.   New->Trans = Copy_Transform(((POLY *)Object)->Trans);
  1480.  
  1481.   for (i = 0; i < term_counts(New->Order); i++)
  1482.   {
  1483.     New->Coeffs[i] = ((POLY *)Object)->Coeffs[i];
  1484.   }
  1485.  
  1486.   return (New);
  1487. }
  1488.  
  1489.  
  1490.  
  1491. /*****************************************************************************
  1492. *
  1493. * FUNCTION
  1494. *
  1495. *   Destroy_Poly
  1496. *
  1497. * INPUT
  1498. *   
  1499. * OUTPUT
  1500. *   
  1501. * RETURNS
  1502. *   
  1503. * AUTHOR
  1504. *
  1505. *   Alexander Enzmann
  1506. *   
  1507. * DESCRIPTION
  1508. *
  1509. *   -
  1510. *
  1511. * CHANGES
  1512. *
  1513. *   -
  1514. *
  1515. ******************************************************************************/
  1516.  
  1517. static void Destroy_Poly(Object)
  1518. OBJECT *Object;
  1519. {
  1520.   Destroy_Transform (((POLY *)Object)->Trans);
  1521.  
  1522.   POV_FREE (((POLY *)Object)->Coeffs);
  1523.  
  1524.   POV_FREE (Object);
  1525. }
  1526.  
  1527.  
  1528.  
  1529. /*****************************************************************************
  1530. *
  1531. * FUNCTION
  1532. *
  1533. *   Compute_Poly_BBox
  1534. *
  1535. * INPUT
  1536. *
  1537. *   Poly - Poly
  1538. *   
  1539. * OUTPUT
  1540. *
  1541. *   Poly
  1542. *
  1543. * RETURNS
  1544. *   
  1545. * AUTHOR
  1546. *
  1547. *   Dieter Bayer
  1548. *   
  1549. * DESCRIPTION
  1550. *
  1551. *   Calculate the bounding box of a poly.
  1552. *
  1553. * CHANGES
  1554. *
  1555. *   Aug 1994 : Creation.
  1556. *
  1557. ******************************************************************************/
  1558.  
  1559. void Compute_Poly_BBox(Poly)
  1560. POLY *Poly;
  1561. {
  1562.   Make_BBox(Poly->BBox, -BOUND_HUGE/2, -BOUND_HUGE/2, -BOUND_HUGE/2, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
  1563.  
  1564.   if (Poly->Clip != NULL)
  1565.   {
  1566.     Poly->BBox = Poly->Clip->BBox;
  1567.   }
  1568. }
  1569.  
  1570.